function [CE,F_star,B_star,Pu_star,Pd_star,F_star_0,B_star_0,Pu_star_0,Pd_star_0] = ...
    Get_CE_TDRILA_optim_fun_parfor(Floor_control,F_choices,Buffer_control,B_choices,Pu_choices,Pd_choices,Pu_Pd_same,u,RILA_fee_rate,c_bar,T,g_I,R_S_GHQ,R_S_GHQ_RN,weight,r,tau,A_grid,I_0)

%% Set model specifications & parameters
if Pu_Pd_same
    Pd_choices = 1;     % we will separately set Pd = Pu during the optimization.
end
F_num = length(F_choices);  B_num = length(B_choices);  Pu_num = length(Pu_choices);    Pd_num = length(Pd_choices);


%% State space
A_num = length(A_grid); A_grid_max = max(A_grid);   % grid of potential QRP account values A_t.


%% Matrizes to store results from dynamic optimization
EU_star = u(0.001) * ones(A_num,T);           % to store expected utility of being in current state (and time)
F_star = zeros(A_num,T-1);          % to store optimal Floor level
B_star = zeros(A_num,T-1);          % to store optimal Buffer level
Pu_star = zeros(A_num,T-1);         % to store optimal Upside Participation Rate
Pd_star = zeros(A_num,T-1);         % to store optimal Downside Participation Rate


%% Terminal utility
EU_star(:,T) = u(max(A_grid,0.01));


%% Recursion
for t = (T-1):-1:1
    EU_tplus1 = EU_star(:,t+1);     % value function at end of period.
    I_t = I_0 * (1+g_I)^t;  % PH's time-t contribution to QRP (no contribution at retirement time T).
    parfor a = 1:A_num
        A_t = A_grid(a);

        % Analyze all potential Floor-based RILA contracts:
        EU_F = u(0.001) * ones(F_num,Pu_num,Pd_num);
        if Floor_control
            for f = 1:F_num
                F = F_choices(f);
                for pu = 1:Pu_num
                    Pu = Pu_choices(pu);
                    for pd = 1:Pd_num
                        if Pu_Pd_same
                            Pd = Pu;
                        else
                            Pd = Pd_choices(pd);
                        end

                        Cap = FindCapRate_GHQ("Floor",F,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                        R_A = max(Pd*R_S_GHQ , -F) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);     % GHQ nodes for return credited to QRP account.
                        A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max );                % GHQ nodes for QRP account value at year-end.
                        U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                      % GHQ nodes for utility at year-end.
                        EU_F(f,pu,pd) = weight' * U_tplus1;                                 % expected utility of account values at end of period.
                    end
                end
            end
        end

        % Analyze all potential Buffer-based RILA contracts:
        EU_B = u(0.001) * ones(B_num,Pu_num,Pd_num);
        if Buffer_control
            EU_B = zeros(B_num,Pu_num,Pd_num);
            for b = 1:B_num
                B = B_choices(b);
                for pu = 1:Pu_num
                    Pu = Pu_choices(pu);
                    for pd = 1:Pd_num
                        if Pu_Pd_same
                            Pd = Pu;
                        else
                            Pd = Pd_choices(pd);
                        end

                        Cap = FindCapRate_GHQ("Buffer",B,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                        R_A = min(Pd*R_S_GHQ+B , 0) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);     % GHQ nodes for return credited to QRP account.
                        A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max );                % GHQ nodes for QRP account value at year-end.
                        U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                      % GHQ nodes for utility at year-end.
                        EU_B(b,pu,pd) = weight' * U_tplus1;                                 % expected utility of account values at end of period.
                    end
                end
            end
        end


        % find combination of (F,B,Pu,Pd) that maximizes expected utility in current state (& time): 
        EU_max_F = max(max(max(EU_F)));
        EU_max_B = max(max(max(EU_B)));
        if EU_max_F > EU_max_B  % Floor preferable to Buffer
            [ff,loc] = find(EU_F == EU_max_F);      % Note: 'loc' will provide the combined location of the largest entry, combining the second and third dimension of the 'EU_B' matrix.
            pdf = floor((loc(1)-1)/Pu_num) + 1;
            puf = loc(1) - (pdf-1)*Pu_num;
            F_star(a,t) = F_choices(ff(1));
            B_star(a,t) = 0;
            Pu_star(a,t) = Pu_choices(puf);
            Pd_star(a,t) = Pd_choices(pdf);
            EU_star(a,t) = EU_max_F;
        else                    % Buffer preferable to Floor
            [bb,loc] = find(EU_B == EU_max_B)      % Note: 'loc' will provide the combined location of the largest entry, combining the second and third dimension of the 'EU_B' matrix.
            pdb = floor((loc(1)-1)/Pu_num) + 1
            pub = loc(1) - (pdb-1)*Pu_num
            F_star(a,t) = 1;
            B_star(a,t) = B_choices(bb(1));
            Pu_star(a,t) = Pu_choices(pub);
            Pd_star(a,t) = Pd_choices(pdb);
            EU_star(a,t) = EU_max_B;
        end

    end
    
end


%% Time 0
EU_tplus1 = EU_star(:,1);     % value function at end of period.
I_t = I_0;      % PH's time-t contribution to QRP (no contribution at retirement time T).
A_t = 0;

% Analyze all potential Floor-based RILA contracts:
EU_F = u(0.001) * ones(F_num,Pu_num,Pd_num);
if Floor_control
    for f = 1:F_num
        F = F_choices(f);
        for pu = 1:Pu_num
            Pu = Pu_choices(pu);
            for pd = 1:Pd_num
                if Pu_Pd_same
                    Pd = Pu;
                else
                    Pd = Pd_choices(pd);
                end

                Cap = FindCapRate_GHQ("Floor",F,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                R_A = max(Pd*R_S_GHQ , -F) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);     % GHQ nodes for return credited to QRP account.
                A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max );                % GHQ nodes for QRP account value at year-end.
                U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                      % GHQ nodes for utility at year-end.
                EU_F(f,pu,pd) = weight' * U_tplus1;                                 % expected utility of account values at end of period.
            end
        end
    end
end

% Analyze all potential Buffer-based RILA contracts:
EU_B = u(0.001) * ones(B_num,Pu_num,Pd_num);
if Buffer_control
    for b = 1:B_num
        B = B_choices(b);
        for pu = 1:Pu_num
            Pu = Pu_choices(pu);
            for pd = 1:Pd_num
                if Pu_Pd_same
                    Pd = Pu;
                else
                    Pd = Pd_choices(pd);
                end

                Cap = FindCapRate_GHQ("Buffer",B,Pu,Pd,r,tau,RILA_fee_rate,c_bar,R_S_GHQ_RN,weight);
                R_A = min(Pd*R_S_GHQ+B , 0) .* (R_S_GHQ<0) + min(Pu*R_S_GHQ,Cap) .* (R_S_GHQ>0);    % GHQ nodes for return credited to QRP account.
                A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-RILA_fee_rate) , A_grid_max );                % GHQ nodes for QRP account value at year-end.
                U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);                      % GHQ nodes for utility at year-end.
                EU_B(b,pu,pd) = weight' * U_tplus1;                                 % expected utility of account values at end of period.
            end
        end
    end
end

% find combination of (F,B,Pu,Pd) that maximizes expected utility in current state (& time):
EU_max_F = max(max(max(EU_F)));
EU_max_B = max(max(max(EU_B)));
if EU_max_F > EU_max_B  % Floor preferable to Buffer
    [ff,loc] = find(EU_F == EU_max_F);      % Note: 'loc' will provide the combined location of the largest entry, combining the second and third dimension of the 'EU_B' matrix.
    pdf = floor((loc(1)-1)/Pu_num) + 1;
    puf = loc(1) - (pdf-1)*Pu_num;
    F_star_0 = F_choices(ff(1));
    B_star_0 = 0;
    Pu_star_0 = Pu_choices(puf);
    Pd_star_0 = Pd_choices(pdf);
    EU_star_0 = EU_max_F;
else                    % Buffer preferable to Floor
    [bb,loc] = find(EU_B == EU_max_B);      % Note: 'loc' will provide the combined location of the largest entry, combining the second and third dimension of the 'EU_B' matrix.
    pdb = floor((loc(1)-1)/Pu_num) + 1;
    pub = loc(1) - (pdb-1)*Pu_num;
    F_star_0 = 1;
    B_star_0 = B_choices(bb(1));
    Pu_star_0 = Pu_choices(pub);
    Pd_star_0 = Pd_choices(pdb);
    EU_star_0 = EU_max_B;
end


%% Convert utility to Certainty Equivalent (CE)
options = optimset('Display','none');
CE = fsolve(@(ce)  1e+3 * (u(ce)/EU_star_0-1), 200,options);    % certainty equivalent to potential values of A_T.